Final 


i 


d/l-J3-gc </y 

r 

/!/ 7 3- S3 *?*/ 


Report for Contract NAS 9-11964 


CASE F' i - E 
COPY 


CALCULATION OF FREE-FALL TRAJECTORIES 
RASED ON NUMERICAL OPTIMIZATION TECHNIQUES 


Prepared for 

National Aeronautics and Space Administration 
Manned Spacecraft Center 
Houston, Texas 


. by 

Texas Center for Research 
Applied Mathematics and Mechanics 
3100 Perry Lane 


Austin, Texas 


78731 


INTRODUCTION 


The purpose of this study has been the development of a 
means of computing free-fall (non- thrusting) trajectories from 
one specified point in the solar system to another specified 
point in the solar system in a given amount of time. The prob- 
lem, therefore, is the determination of the initial velocity 
which--in combination with the initial position and time--can 
be numerically integrated forward in time to the specified 
final time to satisfy the final position requirement to within 
some given tolerance. 

As stated above, the problem is that of solving a two- 
point boundary value problem for which the initial slope is un- 
known. Two standard methods of attack exist for solving two- 
point boundary value problems. 

The first method is known as the initial value or shooting 
method. All unknown quantities at the initial time are guessed 
and the non-linear differential equations of motion are inte- 
grated numerically to the final time. If the final miss dis- 
tance is unsatisfactory, some corrective scheme is applied to 
determine a new value for the guessed quantities and the proc- 
ess is repeated. A disadvantage of this method is that the 
initial values which were guessed must be "fairly close" to the 
correct values or the shooting method may not converge. "Fairly 
close" is a property of the individual problem for non-linear 
problems and it is often most of the work to determine the ini- 
tial guesses which are close enough to converge the initial 
value process. An advantage of this method is that the differ- 
ential equations of motion are satisfied to the tolerance of the 
integration for each iteration and the converged trajectory 
represents a true solution to the integration accuracy. 

The second method of attack for two-point boundary value 
problems is to approximate the non-linear differential equations 
by an appropriate linearized set. The solution method then con- 
sists of solving a series of linear boundary value problems 
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which can be made to satisfy boundary conditions on every iter- 
ate. This series of linear boundary value problems then ap- 
proaches the non-linear problem until the difference between 
the linear solution and the non-linear solution are less than 
some specified tolerance. 

The method of this study uses parts of both boundary value 
problem solution techniques described above. A complete velo- 
city history is guessed such that the corresponding position 
history satisfies the given boundary conditions at the appro- 
priate times. An iterative procedure is then followed until 
the last guessed velocity history and the velocity history 
obtained from integrating the acceleration history agree to some 
specified tolerance everywhere along the trajectory. Conver- 
gence for this method is obtained for fairly poor initial guesses 
of the velocity history and terminal convergence is obtained in 
a quadratic manner. 


DESCRIPTION OF THE METHOD 

Given a two-point boundary value problem for which initial 
position and time and the final position and time are speci- 
fied, determine the initial velocity which can be used in an 
initial value integration scheme to satisfy the final conditions. 

The procedure to be followed is described in detail below. 
First, an initial guess is made of the entire velocity history. 
This guess is made in a manner which satisfies the boundary 
conditions when integrated. Designate this guessed velocity 
history as the control velocity history to differentiate it 
from an integrated velocity history obtained from the integrated 
acceleration. The position is obtained by integrating the con- 
trol velocity history. The position dependent acceleration is 
integrated simultaneously to obtain an integrated velocity his- 
tory. At the final time, two velocity histories are known: 
the control velocity history which satisfies the boundary con- 
ditions, and the integrated velocity history which satisfies the 



3 


non-linear differentia] equations of motion. An iteration 
scheme is now used to drive the control velocity history toward 
satisfaction of the non-linear differential equations while the 
integrated velocity is driven toward satisfaction of the boun- 
dary conditions. The method is considered to be converged when 
the control velocity history and the integrated velocity history 
agree to some specified tolerance everywhere along the path. 

Therefore, a two-part correction scheme is applied to the 
old control velocity history to yield a new control velocity 
history. The first correction assures satisfaction of the boun- 
dary conditions with the new control velocity history and the 
second correction augments convergence to the integrated velo- 
city history from the control velocity history. 


DEVELOPMENT OF THE EQUATIONS 

The formulation of the two-point boundary value problem 
with given initial and final positions and time may be stated in 
the following manner: 


Determine 

V 

ftp 

= v . 
1 

for 




• 

x = 

v(t) 




• 

V = 

F(x,t) 




t i » 

x^ Given 




t £ , 

Xf Given 

where 

X 

is a 

3 x 

1 position matrix 


V 

is a 

3 x 

1 velocity matrix 


F 

is a 

3 x 

1 force matrix 


t 

is a 

scalar 


Guess a velocity history which satisfies the prescribed 
boundary conditions and designate this as u(t) . Also, desig- 
nate the corresponding integrated velocity history as A(t). 
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Now x = u(t) 

X = F(x,t) 

t . , x . Given 
1 1 

t^, x ^ Given 

The last set of equations will be solved on each iteration 
until the new guessed velocity history and the new integrated 
velocity history agree to some specified tolerance. 

The integrated velocity history is to be changed by 6A(t) 
to satisfy the boundary conditions on the next iterate while 
the control velocity history is changed by fiu(t) to agree with 
the integrated velocity history on the next iteration. 

The following relations are desired on the H + 1 st iter- 
ation. 

U £ + l ( ' t ^ = A £ + l^ t ^ 

u *V t3 = Vo * 

Vi Ct) ■ Vo ♦ «vo 

<Su £ (t) = fiA £ (t) + A £ (t) - u £ (t) 

or 5u = fiA + P(A - u) 

The multiplier P in the above expression is a weighting 
matrix which can be preset to some constant value or included 
in the computation. Since P represents the partial derivative 
matrix of the control velocity to the integrated velocity, the 
secant method could be applied to each point to accelerate con- 
vergence. For this study, a constant scalar P where P <_ 1 
was used. 

Linearizing the differential equations leads to the follow- 
ing linear perturbation equations. 

fix = fiu = 5A + P(A - u) 
fiA = F fix 

X 



Assuming that fix and 6 A can be written as linear func 
tions of their initial conditions (with time - dependent coeffi- 
cients) results in the following relations when the fact that 
6x. = 0 is taken into account. 

fix = A (t) + M (t) 

6A = B (t) 5A i + N (t) 


where 

A 

is a 

3x3 

matrix with A. = 0 
i 


M 

is a 

3x1 

matrix with M. = 0 
i 


B 

is a 

3x3 

matrix with B^ = 0 


N 

is a 

3x1 

matrix with N- = 0 
i 

Now, the 

last relations 

are differentiated to yield 



fix = 

A fiX . 

l 

+ M 



fiA = 

B fiX i 

+ N 



fi x = 

6u = 

fiX + P(X - u) 



fiA = 

BfiX i 

+ N + P(X - u) 



6X = 

F fix 

X 




fiX = 

F x AfiX 

. + F M 

IX 

Comparing 

coefficients , 

the following differential equa 


tions for the linear perturbation mapping result. 

A = B 

B = F A 
x 

M = N + p(x - U ) 

N = F M 

X 

These equations are integrated along with the non-linear 
equations. At the final time: 
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<5Xf = x f - x(t f ) 

6x £ = A r 6A- + Mr 

<SX ± = A^ 1 (x f - x(t £ ) - M f ) 

This provides the new value for the initial velocity. 

S A = B6A- + N 

l 

6u = 6A + P ( A - u) 

<5u = B6A. + N + P ( A - u) 

6u . = 6 A . 

l l 

This is the procedure to follow in making changes in the 
old control velocity history to begin a new iteration. The 
iteration procedure is stopped when the maximum value of Su 
anywhere along the path is less than the specified tolerance. 

NUMERICAL INTEGRATION AND INTERPOLATION 

In order to apply the iteration scheme in the preceeding 
section, it would be advantageous to employ a constant stepsize 
numerical integration procedure so that the corrections to the 
control velocity history would always occur at the tabulated 
points of the control velocity history. The terminal conver- 
gence properties and the integration accuracy are limited in a 
constant stepsize integrator- -particularly if any singularities 
are present in the interval of interest. 

Due to these limitations, a variable step Runge-Kutta- 
Fehlberg integrator with truncation error control was chosen 
for the integration package. In order to use this integrator, 
an interpolation scheme for the control velocity history is 
necessary to provide values of the control velocity between the 
tabulated points on each iteration. The stepsize pattern 
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changes from iteration to iteration so the tabulated control 
velocity points also change from iteration to iteration. For 
the interpolator, a cubic spline was chosen. The cubic spline 
fits a third order polynomial in time between tabulated points 
and provides continuous function values and continuous first 
and second derivative values throughout the entire control ve- 
locity history. The cubic spline was chosen for its stability, 
availability, and applicability. 

With the cubic spline fitting a third order polynomial in 
each stepsize interval, the natural choice of the integrator 
order was determined to be the RKF 3(4) integrator which inte- 
grates a third order polynomial to the specified integration 
tolerance. The integrator, therefore, chooses a stepsize for 
a third order polynomial on integration and the spline fits a 
third order polynomial of interpolation to this interval. 
Although the combination of RKF 3(4) integrator and cubic spline 
should be optimum for accuracy, the combination is not the opti- 
mum for minimizing storage requirements. A higher order inte- 
grator takes fewer steps for the same integration and, there- 
fore, requires less storage for the tabulated velocity history. 
The spline fit to this larger interval is not as accurate as in 
the preceeding combination. A compromise is, therefore, indi- 
cated- -the choice of integrator to couple to the cubic spline 
being governed by the accuracy and storage requirements. An 
optimum compromise between accuracy and storage is probably the 
RKF 4(5) integrator coupled to the cubic spline. 

For the initial value process which follows the conver- 
gence of the control velocity iteration, a high order integrator 
[RKF 7(8)] with a fine tolerance is used to refine the trajec- 
tory. The high order integrator is used in this part for high 
accuracy with few steps. Since this is an initial value proc- 
ess, no interpolation is necessary and the differential equa- 
tions are satisfied to the integration accuracy. 
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APPLICATIONS 

The initial test of the control velocity iteration scheme 
was chosen to be the generation of Apollo-type earth-moon tra- 
jectories. The initial trajectories have been generated for 
both coplanar and non-coplanar cases. The coplanar cases util- 
ize initial and final points which lie in the moon's orbital 
plane and the non-coplanar cases have final points which do not 
lie in this plane. The initial point was chosen for a 100 mile 
departure altitude from earth and a 50 mile arrival altitude at 
the moon. 

Beginning with the coplanar cases, a control velocity his- 
tory was chosen which kept the spacecraft in the lunar orbital 
plane at all times. The initial guess missed the desired final 
point by over 2000 miles. Convergence of the control velocity 
iteration scheme required 9 iterations for a total of less than 
20 seconds computation time on the CDC 6600. Using this ini- 
tial velocity in a standard perturbation scheme (initial value 
method), convergence to a position miss of less than 1 inch re- 
quired 3 iterations for a total of 7 additional seconds. 

Choosing a non-coplanar initial velocity history guess with 
the same initial and final points as the case above, the control 
velocity iteration scheme converged back to the coplanar veloc- 
ity history. Again, the perturbation scheme converged in 3 
iterations. Run times for both sections were comparable to the 
run times of the first case. 

Moving the final point 24 miles out of plane caused the 
third case to converge to a non-coplanar velocity distribution 
f rom the initial guess of the first case. Run times were again 
compatible with the first case. 

The control velocity iteration scheme was next applied to 
the generation of interplanetary trajectories. Initially, the 
planets were placed in circular orbits about a fixed sun. 

Several trajectories of the earth-mars, earth-venus, earth- 
mercury types have been generated for short flight times 
on the order of 50 to 100 days. The short flight time 
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allows the use of a constant velocity for the initial guess at 
the velocity history. The trajectories have been generated for 
an initial point which is at a 300 mile altitude above the earth 
and a final point which is 300 miles altitude at the target 
planet. During the computation of these trajectories, the sun 
and the five inner planets were active gravitational sources. 
Computation time and accuracy for these trajectories were com- 
parable to those of the earth-moon trajectories. For example, 
a 90-day earth-venus trajectory from 300 mile altitude to 300 
mile altitude converged the control velocity iteration scheme in 
5 iterations for a computation time of less than 24 seconds. 

The perturbation scheme again converged in 3 iterations. 

A second planetary ephemeris which used 3 - dimens ional el- 
liptic orbits with constant elements was the next improvement. 

An attempt was then made to generate an earth-saturn trajectory 
with an intermediate jupiter flyby. The control velocity iter- 
ation scheme converged to a trajectory from the initial point 
to the final point but the flyby at jupiter was not present. 

The initial corrections to the velocity history move the trajec- 
tory path away from the flyby local extremum to the strong non- 
flyby extremum. A change in the flight time might bring this 
trajectory back within the influence of jupiter to permit a 
flyby . 

Several other short flight time trajectories between two 
planets were converged and it was observed that the 3-dimen- 
sional trajectories were more difficult to converge than the 
2 - dimens ional trajectories. 

The third planetary ephemeris to be incorporated was the 
J.P.L. Analytical ephemeris based upon polynomial approxima- 
tions to the orbital elements of the planets. This ephemeris 
was chosen as the final one and a specific trajectory was 
chosen for computation. The specific mission was an earth 
to venus trajectory with a heliocentric transfer angle greater 
than 200°. This mission was chosen for its importance in the 
radioactive waste disposal project. 
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The trajectory end points were chosen such that the initial 
point was at a 100 mile altitude over a specific point on the 
earth and the final point was chosen to be at a 100 mile alti- 
tude over a specific point on venus. The earth launch point 
was chosen for compatibility with existing launch facilities 
and the venus final point was chosen from preliminary patched 
conic studies . 

The initial velocity guess would ideally be the velocity 
distribution from the patched conic trajectory which would con- 
sist of hyperbolic planetocentric segments matched in position 
and velocity to elliptical heliocentric segments. Unfortunately, 
this system was not available. The result was that Lambert’s 
theorem was used to generate a heliocentric elliptical segment 
connecting the initial and final points. If the planetary grav- 
itational fields are introduced full strength, the trajectory 
integration requires more integration steps than the allotted 
storage permits. 

The reason for this was that the elliptic section which was 
the solution of Lambert's theorem passed through the earth near 
to its center. The variable step integrator, therefore, needed 
more and smaller steps for an accurate integration. To circum- 
vent this difficulty, the planet's gravitational strength was 
gradually increased from zero to 1001 in steps. Lambert's 
theorem provided the solution for 01, then a II problem was con- 
verged, then a 101 problem and finally a 1001 problem. This 
allows the control velocity iteration scheme to shape the tra- 
jectory gradually. More time is required for this procedure 
but final convergence is obtained without exceeding the allotted 
storage. 

Beginning with the Lambert's theorem velocity history 
guess, 4 iterations were required for convergence of the II 
problem, 6 iterations were required for convergence of the 101 
problem, and 10 iterations were required for convergence of the 
1001 problem. The total run time was 179.59 seconds. 

A better initial velocity guess history would do away with 
the relaxation of the planetary gravitational fields and also 
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augment convergence. This would greatly reduce the total re- 
quired run time. 

The initial velocity was taken from the above converged 
run and used in a standard initial value perturbation scheme. 

The integrated terminal miss was reduced somewhat but quad- 
ratic convergence was not obtained. The initial integration of 
this scheme usually shows a greater terminal miss than the con- 
verged velocity guess scheme since the differential equations 
of motion are obeyed exactly in the standard perturbation scheme. 
A reduce- the-norm iteration procedure was employed with the 
standard perturbation scheme which would not accept an iterate 
unless it reduced the norm of the miss distance from the last 
accepted iterate. 

To generate a flyby trajectory, two problems of the pro- 
ceeding type were solved. In the velocity guess iteration 
scheme, the joining point was fixed at venus and the two trajec- 
tories (earth- to -venus and venus-to-terminal-point) were matched 
in position (to integration accuracy) but not in velocity at the 
joining point. An iteration scheme was then employed in the 
standard perturbation scheme which was to iterate out the veloc- 
ity of each segment to vary. This process was successful in 
decreasing the velocity mismatch to 30-40 fps but integration 
accuracy was such that no further reduction could be obtained. 


CONCLUSIONS AND RECOMMENDATIONS 

The control velocity iteration is a very powerful tool for 
solving two-point boundary value problems. Convergence can be 
obtained to the accuracy of the interpolation scheme even for a 
somewhat poor initial guess. The method is particularly appli- 
cable to the computation of velocity histories for interplanetary 
trajectories where the initial velocity guess might not con- 
verge the standard perturbation scheme. 

The generation of planetary flyby trajectories by the veloc- 
ity guess iteration will require the inclusion of a constraint 
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at the flyby planet in order to assure that the trajectory con- 
verges to the flyby which is a weak local extremum instead of 
the non- flyby which is a much stronger extremum. Without this 
constraint, the early stages of the velocity iteration scheme 
move the trajectory away from the flyby trajectory to a non- 
flyby trajectory. 

The accuracy of the integration should be improved in the 
neighborhood of the planets if the integration scheme is switched 
from heliocentric ecliptic to planetocentric ecliptic when inside 
the planet's sphere of influence. More accuracy of the integra- 
tor and, therefore, of the interpolator should improve the con- 
vergence characteristics of the velocity guess iteration. The 
accuracy should also be reflected in the initial value perturba- 
tion scheme by providing much more accurate information for com- 
puting changes in the initial velocity. 



